pacmanを使うと、install.packages()やlibrary()を使わずにパッケージをまとめて読み込むことができる。
まずpackmanをインストールする。
install.packages("pacman")今回は、readxlとtidyverseとplotlyとskimrを使う。以下のように書けば、まとめて読み込むことができる。
pacman::p_load(readxl, tidyverse, plotly, skimr)作業ディレクトリを設定する。プロジェクト機能を活用できるのであれば、そちらを活用する。
setwd(自分のディレクトリ) # ディレクトリの設定
getwd # ディレクトリの確認前回同様、今回も test_scores.xlsx のデータを使う
df <- read_excel("test_scores.xlsx") # エクセル読み込み
head(df)## # A tibble: 6 × 5
## クラス 名前 数学 英語 国語
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 C 安藤 67 78 45
## 2 A 藤村 55 73 72
## 3 B 大泉 47 55 73
## 4 B 木村 78 71 85
## 5 A 鈴井 51 74 61
## 6 C 望月 79 84 75
今回は、dplyr::rename()によって変数名を英語に変換してから使う。
#dfの名前変更
df <- df %>% dplyr::rename(class = クラス,
name = 名前,
math = 数学,
english = 英語,
japanese = 国語)
head(df)## # A tibble: 6 × 5
## class name math english japanese
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 C 安藤 67 78 45
## 2 A 藤村 55 73 72
## 3 B 大泉 47 55 73
## 4 B 木村 78 71 85
## 5 A 鈴井 51 74 61
## 6 C 望月 79 84 75
グラフの描画で役立つパッケージを紹介する。
{ggplot2}というグラフの描画用パッケージは綺麗なグラフを描くことができる。
ggplotは次のように2つ以上の手順をとり、2つ以上の関数を+でつなげていき、レイヤーを重ねるように描画していく。
ggplot(データフレーム, aes(x, y)))geom_histogram()などgeom_なんたらという関数)※R Studio Cloudで図を描くとき、タイトルや軸タイトルに日本語を使うと文字化けしてしまう。したがって、英語やローマ字を使うか、日本語を使いたい場合はパソコンにRとR Studioをインストールして使用する必要がある。(2019.4.22時点)
# 棒グラフをつくり、gに代入する
g <- ggplot(df, aes(x = class)) + # dfのデータでキャンバスを作る
geom_bar() + # キャンバスに棒グラフgeom_bar()を重ねる
xlab("Class") # X軸タイトルの変更(R Studio Cloudでは日本語文字化けするため)
g # 表示{plotly}パッケージを使うと、簡単なインタラクティブなグラフ(カーソルを合わせると情報が表示されたり、表示情報を変更したりできるグラフ)を描くことができる。
ggplotly()関数を使うと、{ggplot2}で描いたグラフをインタラクティブなものに変換できる。
※R Studio
Cloudで {ggplotly} を使うと、日本語表示は
{ggplot2}
そのものよりも適切に表示されるが、やはり少しおかしい。(2019.4.22時点)
plotly::ggplotly(g) # ggplot2のグラフをplotlyのグラフへ変換データの分布を要約するための指標(要約統計量)の代表的なものが平均(mean)と分散(variance)である。
\[ \bar{x} = \frac{1}{n}\sum_{i=1}^n x_i \]
\[ \sigma^2=\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^2 \]
ただし、Rのvar()関数の算出する値は、母集団の分散(母分散)の不偏推定量としての分散すなわち不偏分散である。
\[ s^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2 \]
# 不偏分散
var(df$english)## [1] 333.1481
分散の単位は元の単位の2乗であるため、それを元の単位に戻すために平方根をとったものが標準偏差(standard deviation)である。
# 不偏標準偏差
sd(df$english)## [1] 18.25234
summary()関数は、データフレームの各行に関する代表値などの情報を取得することができる関数である。
summary(df)## class name math english
## Length:40 Length:40 Min. : 42.00 Min. :42.00
## Class :character Class :character 1st Qu.: 51.75 1st Qu.:50.75
## Mode :character Mode :character Median : 65.00 Median :65.50
## Mean : 66.58 Mean :68.33
## 3rd Qu.: 80.50 3rd Qu.:82.25
## Max. :100.00 Max. :99.00
## japanese
## Min. : 40.00
## 1st Qu.: 51.00
## Median : 67.50
## Mean : 66.88
## 3rd Qu.: 79.50
## Max. :100.00
classとnameは文字列(character)型なので、Length(行数、標本数)が40であることくらいしか表示されない。
english、japanese、mathは整数(integer)型なので、代表値などの要約統計量(summary
statistics)が表示されている。
それぞれの統計量の意味は
Min.が最小値1st Qu.が第一四分位数(1st
quantile):「データを小さいものから順に並べた時に標本数の4分の1の位置(下位25%)にくる値」Medianが中央値(median):「データを小さい順に並べた時に標本数の2分の1の位置にくる値」Meanが平均値(mean):「データの総計をデータの個数で割ったもの」\(\sum_{i=1}^n x_i/n\)3rd Qu.が第三四分位数(3rd
quantile):「データを小さい順に並べた時に標本数の4分の3の位置(75%)にくる値」Max.が最大値である。これらのうちMean以外の5つの統計量をまとめて五数要約という。
平均値は外れ値(極端に高かったり低かったりする値)の影響を受けやすいため、これらの統計量を見ることで、中央値と平均値が大きく異なるデータは外れ値があったりデータの分布に歪みがあることが推測できるわけである。
skimrのskim()を使うと、小さなヒストグラムつきの要約統計量の表を出力できる。欠損値の数なども表示してくれる。
skimr::skim(df)| Name | df |
| Number of rows | 40 |
| Number of columns | 5 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| class | 0 | 1 | 1 | 1 | 0 | 3 | 0 |
| name | 0 | 1 | 1 | 3 | 0 | 38 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| math | 0 | 1 | 66.58 | 18.13 | 42 | 51.75 | 65.0 | 80.50 | 100 | ▇▆▅▃▅ |
| english | 0 | 1 | 68.33 | 18.25 | 42 | 50.75 | 65.5 | 82.25 | 99 | ▇▅▂▆▅ |
| japanese | 0 | 1 | 66.88 | 17.17 | 40 | 51.00 | 67.5 | 79.50 | 100 | ▇▃▇▅▃ |
相関係数(correlation coefficient)は2つの量的変数の直線的な関係の強さを測るもので、次の式で定義される
\[ \begin{align} x\text{と}y\text{の相関係数}r &= \frac{x\text{と}y\text{の共分散}}{x\text{の標準偏差}\times y\text{の標準偏差}} \\ &= \frac{\frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})(y_i - \bar{y})} {\sqrt{\frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2} \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - \bar{y})^2}} \end{align} \]
Rではcor()で計算できる。相関係数は-1から+1の範囲の値をとり、絶対値で1に近いほど関係が強いことを示す。
cor(x = df$english, y = df$math)## [1] 0.160109
相関係数が入った行列を相関行列という
cor(df[,3:5])## math english japanese
## math 1.000000 0.16010896 0.22593799
## english 0.160109 1.00000000 0.04096556
## japanese 0.225938 0.04096556 1.00000000
ヒストグラムは量的変数の分布を確認するには極めて有用なグラフである。
# ヒストグラム
ggplot(df, aes(x = english)) +
geom_histogram(bins = 10) + #指定しない場合はbins binsは棒の数。=30になる
xlab("English") #X軸タイトルを変更(R Studio cloudは日本語文字化けするため)R本体の機能でヒストグラムを描く場合はhist()を使う
# ヒストグラム
hist(x = df[["math"]])# ggplotの棒グラフ
g <- ggplot(df, aes(x = class)) +
geom_bar() +
xlab("Class") #X軸タイトルを変更
#グラフ表示
g#plolyで表示
plotly::ggplotly(g) 2変数でグラフを描く場合、キャンバスにx(グラフ横軸)とy(グラフ縦軸)の両方を指定する。
そして散布図を描く場合はgeom_point()を+でつなげればよい
# ggplot
g <- ggplot(df, aes(x = english, y = japanese)) +
geom_point() +
xlab("English") + #X軸タイトルを変更
ylab("Japanese") #y軸タイトルを変更
# 散布図表示
g# plotlyで表示
plotly::ggplotly(g) 回帰直線を加えるには、stat_smooth()を使う。
ggplot(df, aes(x = english, y = japanese)) +
geom_point() +
stat_smooth(method = "lm", se =FALSE, colour = "blue") + # lm:最小二乗法、信頼区間は表示させない(FALSE)、色は青(blue)
xlab("English") + #X軸タイトルを変更
ylab("Japanese") #y軸タイトルを変更折れ線グラフはgeom_line()で描くことができる。
Excelで言うところの「マーカー付き折れ線グラフ」にしたい場合はgeom_point()を重ねればよい。
ggplot(df, aes(x = english, y = japanese)) +
geom_line()+ # 折れ線
geom_point() + # マーカー
xlab("English") +
ylab("Japanese")上では説明の簡単のため散布図と同じデータを使ったが、折れ線グラフは本来は時系列データに使うのが望ましい。
そこで、以下ではRにあらかじめ収録されている練習用のデータセットのひとつであるlongleyを使った例をのせる。
このデータは1947~62年のアメリカのGNPや雇用者数などが収録されている(詳細はhelpを参照してほしい)。
head(longley)## GNP.deflator GNP Unemployed Armed.Forces Population Year Employed
## 1947 83.0 234.289 235.6 159.0 107.608 1947 60.323
## 1948 88.5 259.426 232.5 145.6 108.632 1948 61.122
## 1949 88.2 258.054 368.2 161.6 109.773 1949 60.171
## 1950 89.5 284.599 335.1 165.0 110.929 1950 61.187
## 1951 96.2 328.975 209.9 309.9 112.075 1951 63.221
## 1952 98.1 346.999 193.2 359.4 113.270 1952 63.639
# ggplot
g <- ggplot(longley, aes(x = Year, y = GNP)) +
geom_line()+ # 折れ線
geom_point()+ # マーカー
labs(title = "1947~62年のアメリカのGNP")
# グラフ表示
g# plotlyで表示
ggplotly(g)箱ひげ図は五数要約の統計量(最小値、第一四分位数、中央値、第三四分位数、最大値)をプロットすることで分布の概形を表示する図である。
g <- ggplot(df, aes(x = class, y = japanese)) +
geom_boxplot()+ # 箱ひげ図
labs(title = "クラスごとの国語の点数の分布") +
xlab("Class") +
ylab("Japanese")
ggplotly(g)バイオリンプロットは、「カーネル密度推定」という手法で推定されたデータ分布の曲線を左右対称に描いたもので、異なるデータの分布を視覚的に比較できる。デフォルトでは最小値から最大値の範囲でプロットされ、geom_violin(trim = FALSE)にすると曲線の端はカットされない。
g <- ggplot(df, aes(x = class, y = japanese)) +
geom_violin() + # バイオリンプロット
labs(title = "クラスごとの国語の点数の分布") +
xlab("Class") +
ylab("Japanese")
ggplotly(g)g <- ggplot(df, aes(x = class, y = japanese)) +
geom_violin() + # バイオリンプロット。これを先に描く。
geom_boxplot(width = 0.3) + # 箱ひげ図を上書き。幅を0.3に指定
labs(title = "クラスごとの国語の点数の分布") +
xlab("Class") +
ylab("Japanese")
g #ggplotly(g)だと箱ひげ図が表示されないggplot()関数の中に入れるaes()関数のfillやcolorといった引数に質的変数を指定することで、色の塗り分けによる層別化ができる。
これにより3つ目の変数の情報を表現できる。
ggplot(df, aes(x = english, y = japanese, fill = class, color = class)) +
geom_point()+
labs(title = "英語と国語の得点の散布図をクラスごとに塗り分けたもの") +
xlab("English") +
ylab("Japanese")散布図のサークルの大きさを、第三の変数に合わせて変更することもできる。
# ggplot
g <- ggplot(df, aes(x = english, y = japanese, size = math)) +
geom_point(alpha = 0.5) + # alpha:バブルの透過度
labs(title = "英語と国語の得点の散布図(バブルサイズは数学の点数を反映)") +
scale_size(range = c(1,10)) + # バブルサイズの範囲を設定
xlab("English") +
ylab("Japanese")
# 散布図表示
g#plotly
ggplotly(g) 図の保存にはggsaveを使用する
ggsave(plot = g, file = "bubble_chart.png") # violin_box.pngとして保存される(必要に応じて保存場所は指定する)データ可視化については、様々な情報やテキストがある。また、グーグルの画像検索から自分が作りたいグラフ・イメージを探し、そのRコードを参考にするのも一つの方法である。
Rob Kabacoff(2020) Data Visualization with R
キーラン・ヒーリー(2021) [『データ分析のためのデータ可視化入門』] (https://www.kspub.co.jp/book/detail/5164044.html) 講談社
ill-identified diary (2021) 『データ分析のためのデータ可視化入門』と最近の R グラフィックスパッケージ事情